Cluster evolution in steady-state two-phase flow in porous media 
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We report numerical studies of the cluster development of two-phase flow in a steady-state en- 
vironment of porous media. This is done by including biperiodic boundary conditions in a two- 
dimensional flow simulator. Initial transients of wetting and non-wetting phases that evolve before 
steady-state has occurred, undergo a cross-over where every initial patterns are broken up. For flow 
dominated by capillary effects with capillary numbers in order of 10~^, we find that around a criti- 
cal saturation of non-wetting fluid the non-wetting clusters of size s have a power-law distribution 
Us ~ s~'^ with the exponent r = 1.92 ± 0.04 for large clusters. This is a lower value than the result 
for ordinary percolation. We also present scaling relation and time evolution of the structure and 
global pressure. 

PACS numbers: 05.40.-a, 07.05. Tp, 47.54.-|-r, 47.55. Mh 



I. INTRODUCTION 

The complex nature of multi-phase flow in porous me- 
dia has acted as motivation for extensive studies in the 
recent years. Different types of fluid displacements in 
porous media play important roles in natural processes, 
many of them of large practical importance such as oil 
recovery, soil mechanics and hydrology. Large quantities 
of petroleum and water resources are hidden in fractured 
and porous stones and that is of vast economic and social 
importance. 

Two-phase displacement in porous media has been 
studied through seminal experimental work Q, |^ |^ , nu- 
merical simulations ^ ^ ,6|| and theoretical work includ- 
ing statistical models and differential equations 0, 
Dynamics of instabilities in immiscible two-phase flow 
controlled by the interplay between viscous and capillary 
forces, give rise to complex pattern formations. 

The traditional way of describing those situations has 
been through either drainage or imbibition. That will 
give rise to transient effects of front propagation such 
as invasion percolation (IP), viscous fingering and stable 
front displacement. The common feature of these effects 
are that they are out of equilibrium, and the IP regime 
will eventually even end with a static regime with im- 
mobile structures due to large capillary barriers. Most 
experiments and simulations in this area deal with sys- 
tems where one and only one fluid is injected into an- 
other until the injected fluid reaches the sink of the sys- 
tem. In natural reservoirs of fluid in porous media the 
situation is, however, dynamic and governed by the inter- 
play and competition between drainage and imbibition. 
This is best described as a steady-state regime inside a 
representative volume element where both drainage and 
imbibition occurs. Pure drainage or pure imbibition are 
pore-scale concepts, and can therefore not alone be scaled 



* Thomas. Ramstad@phys.ntnu.no 
tAlex.Hansen@ntnu.no 



up in a steady-state situation. 

In this paper we examine distribution of cluster for- 
mations in a pore scale network under both steady-state 
and transients. This study contains both quantitative 
and qualitative results of cluster distribution and dy- 
namic evolution of the phases. The model used in our 
simulations is mainly similar to that proposed by Aker 
et al 9] with later generalizations by Knudsen et al [loj 
to include biperiodic boundary conditions. The bound- 
ary conditions in the flow direction lets the fluid that 
leaves the system enter the system again at the inlet of 
the model in a seamless manner. 

This ables us to investigate steady-state flow which in- 
flects the situation deep inside reservoirs. With both an 
invading front of non-wetting fluid and wetting imbibi- 
tion the system will give rise to a blending of structure 
formations. These formations are formed by the interplay 
of different effects like fragmentation of large clusters, 
merging of smaller clusters and a diffusion of fragments. 

This paper is organized in the following way. In Sec- 
tion we describe outlines of the model used for our 
simulations, and in Section IlIII we sketch the initial con- 
ditions used followed by a discussion of the qualitative 
behaviour of initial transients. 

Further, we discuss quantitative features of the simu- 
lations in Section [l VI prior to a scaling analysis. 



II. MODEL 

The main statistical models that reproduce the basic 
domains of porous flow belong to the family of growth 
models that obey the Laplacian equation V^P — Q where 
P is a pressure field and with an interfacial growth rate 
q oc VP. 

The simulations used in this work is based upon a basis 
model of disordered media, that consists of tubes orien- 
tated 45 degrees to the overall flow direction •^,1^. The 
network has a coordination number of 4, and the arrange- 
ment of the tubes is square planar. The tubes consist of 
both the pore volume and the throat volume. The points 
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where the tubes meet are referred to as nodes. Random- 
ness is incorporated in the model by allowing the length 
of the tubes dij to be chosen randomly within 30% of the 
mean length i.e. the lattice constant d, and the radii to 
be chosen from a flat distribution r g (O.ld, O.ld+O.Sdy). 

Distributions of fluid that we report are highly depen- 
dent of saturation of phases and Ca. The capillary num- 
ber 



Ca = 



fJ-Qt 



7S 



(1) 



denotes the ratio between viscous and capillary forces 
where fj, is the largest viscosity rate of the two fluids, 
7 is the surface tension, S denotes the total input sur- 
face of the system and Qtot is the global flow rate in the 
system. The patterns that form during invasion of a non- 
wetting fluid into a defending wetting one are character- 
ized mainly into viscous fingering, capillary fingering or 
stable displacement depending on Ca and the viscosity 
ratio of the defending fii and invading phase /X2 given by 



M = 



(2) 



Our system of tubes is filled with two fiuids, wetting 
and non-wetting respectively, each having a certain frac- 
tion of the total volume denoted non wetting saturation 
and wetting saturation Sv,- The fluids are separated 
through many menisci. The model does not include film 
flow and only one bubble is allowed inside one tube at 
the time. This means that if more than two menisci are 
created within a tube during the same time step, they 
are subsequently collapsed into two menisci. 

Since the fluids are immiscible, the meniscus gives rise 
to a capillary pressure within the tube. The capillary 
pressure is dependent on both the surface tension 7 and 
the radius of the tube given by the Young-Laplace law 



27 , 

Pc = COS f 

r 



(3) 



where 9 is the wetting angle, i.e. the contact angle be- 
tween the wetting phase and the wall of the tube. In our 
model, the tubes are considered to be cylindrical with 
respect to permeability, but hour-glass shaped with re- 
spect to capillary pressure. The relation for the capillary 
pressure reads 



27 



Pc — — (1 — cos(27rx)) 
r 



(4) 



where x is the normalized position of the meniscus run- 
ning from to 1. The local flow rate g of a tube follows 
the Washburn equation for capillary flow , 



k nr 
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{Ap - Pc 



(5) 



where pc is given by Eq. Q and the permeability k = 
This permeability is taken directly from Hagen- 
PoiseuUe flow in cylindrical tubes. We define the effective 
viscosity as 



(6) 



with respect to the position of the non- wetting meniscus 
Xnw The pressure difference between the inlet and outlet 
of a tube is Ap. 

To solve the transport equations we use the property 
that no net volume can be stored in the nodes connecting 
tubes. The net flux through a node is therefore zero and 
this gives a set of linear equations to solve. To find the 
positions of menisci in the tubes we integrate Eq. (O 
forward with a predefined time step, Ai, according to an 
explicit Euler scheme. Inside a tube one meniscus moves 
with a front speed determined from the local flow rate 
q, but when a meniscus reaches the end of a tube it is 
further distributed among the other tubes with ingoing 
flux. Here it is crucial that volume is conserved. Further 
details can be found in ref. 0, . 

The majority of existing models are quasi static and 
therefore apply only when the capillary forces are domi- 
nant. This model, however, accounts for dynamic effects 
and is therefore also capable of handling fast flow or flow 
with little surface tension, hence a wider range of the 
capillary number Ca. 



A. Boundary conditions 

Most experimental setups and simulations are done on 
systems that are out of steady-state. They mainly deal 
with systems where one fluid is injected into another de- 
fending fluid and ends when the invading fluid as reached 
the sink of the model. 

In order to simulate steady-state flow, we use biperiodic 
boundary conditions. This is done by connecting the inlet 
and outlet row, placing the system on a surface with one 
extra dimension. In two dimensions this is easily pictured 
by placing it on a three-dimensional torus as shown in 
Fig. [2 

To make the system develop in time, the global pres- 
sure field AP is applied over the row where the original 
inlet and outlet rows meet. This is done by the use of 
ghost sites jl^ where the global pressure is applied. 

With biperiodicy in the boundary conditions we will 
have two invading processes that compete. One being the 
original invading process as described above with a non- 
wetting phase penetrating a wetting phase, but there will 
also be a wetting phase displacing a non-wetting phase 
since the wetting fluid imbibition from the system will 
enter it at the inlet row. 
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Figure 1: Schematic figure that shows the interpretation the biperiodic boundary conditions mean to our simulations. In order 
to model a situation deep inside a natural reservoir in steady state, we must consider the system size as a small part of a greater 
system. Statistically, it is the same fluid-configurations entering the element as are leaving it. This is illustrated in the figure 
on the left hand side. The way we consider our system is shown on the figure to the right where the two-dimensional system 
is mapped on to a three-dimension torus. 



III. INITIAL TRANSIENTS 

In this section we describe the qualitative behaviour of 
our system as it evolves. As time progresses, the dynamic 
evolution in our model will give rise to different transient 
effects that eventuahy will lead to a steady-state config- 
uration. These transients are dependent on the control 
parameters in our system. 

Typical parameters that are fixed during our simula- 
tions are surface tension 7 and viscosity /i. We keep 
7 — 30.0 mN/m and /inw — ^-w — O.I Pa s which gives 
the viscosity ratio M = 1. In order to set the capillary 
number, we control the overall How rate Qtot- 

In this paper we consider networks of sizes up to 
1024 X 1024 nodes. We wish to analyze the distributions 
of clusters with different Ca and different system sizes. 

At the beginning of the simulations, the network is seg- 
regated into one part of non-wetting fluid and one part 
containing only wetting fluid. At the start of the simula- 
tion the non- wetting fluid will start to invade the wetting 
fluid part. Depending on the capillary number it will ei- 
ther start a viscous fingering behaviour for high Ca or an 
IP process when the capillary forces dominate. The latter 
case is visualized in Fig. [5] where the drainage takes place 
in the upper part of the system. In the early stages pre- 
sented in the figure, a "belt" of non-wetting (black) fluid 
spans the system in horizontal direction. This "belt" is 
an initial configuration chosen by us and will be referred 
to later in this paper. 

On the other hand, wetting fluid will enter the original 
non-wetting region from below pushing the non-wetting 
fluid in the advancing direction. This displacement of 
non-wetting fluid in favour of the wetting happens in 



a piston like manner. There are some trapped regions 
left of non-wetting fluid, but these are smaller than the 
trapped wetting clusters in the other end of the model. 
Since the non- wetting front fingers advance more rapidly 
than the wetting, more compact propagating front |29l| . 
it will sooner or later catch up with the front of the ad- 
vancing wetting fluid. At this point there will be a com- 
petition between the two different advancing fronts. 

When the system reaches steady-state there will be a 
mixing of the two phases. This mixing is dependent both 
of Ca and the saturation of fluids. If the difference in the 
saturation of fluids is large one of the phases will span 
the system entirely, in steady-state, with a large cluster. 
How large this difference needs to be in order to get a 
spanning cluster, depends on the capillary forces and the 
structure of the network i.e. the porosity. If the capillary 
forces dominate over the viscous ones, the single phase 
pressure Ps according to Darcy's law 

Q ^ kA^ 

where L is the system length along the single phase gra- 
dient APg, will be significantly smaller than the total 
pressure gradient. In this case if there is a spanning clus- 
ter of one fluid it will be much easier for the fluid to move 
inside this cluster because of the high capillary pressure 
threshold in non-occupied tubes. As a consequence once 
a spanning cluster evolves, the total pressure will drop 
and one phase stays immobile or trapped. 

Since the total flux is constant and the capillary pres- 
sure will vary with time as menisci configurations change, 
the global pressure will fluctuate. However, when the 
system reaches steady-state the pressure will fluctuate 
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Figure 2: (Color online) Four different stages of a low Ca flow configuration with viscosity match M — 1.0. Ca — 3.2 x 10~^ 
implies a capillary-dominant regime. The different snapshots of the ffow patterns are taken at t = 5000s, t — 10000s, t = 15000s 
and t = 25000s. The system size is 128 x 128 with a Snw ~ 0.5. The largest connected cluster is colored red in the online 
version of the article. The two first snapshot corresponds to the first regime of the pressure in Fig. |3 In the third snapshot 
some of the initial configuration is still left and the pressure is building up according to Fig. |3] The last one is in the steady 
state regime where every initial structure has vanished. 



around a mean value. How long it takes before our sys- 
tem reaches steady-state is a matter of the initial con- 
figuration of the two fluids. As the pressure relaxes 
around a mean value the fractional flow of each fluid 
will settle. The non wetting fractional flow is defined 
as Faw = Qnw/Qtot and the wetting fractional flow is 

It has been shown in previous work m m that i^nw 
depends on both the capillary number and the saturation 
S'nw of the non-wetting fluid. For high values of Ca, the 
curve of i^nw will behave almost linearly with respect to 
S'nw, while in the case of very small Ca we will get regimes 
where one fluid is almost immobile. This feature was also 



found to be independent of system size. 

For steady-state flow we are only able to obtain a to- 
tally immobile regime or a dynamic regime where clus- 
ters rearrange themselves continuously. The first situa- 
tion can be characterized as traditional invasion percola- 
tion where capillary forces dominate. When a percolat- 
ing cluster arises, the structure stays static as long as no 
global parameters are changed. There will be a sudden 
drop in global pressure P as no capillary barriers will 
have to be overcome to sustain the global flow rate. 

In the opposite case, clusters of fluid move, fragment 
and coalesce in an equilibrium process with a global pres- 
sure fluctuating around a mean value. This is true for 
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the viscous case where clusters are not trapped inside a 
growing percolating cluster. However, we report in this 
paper that this is also the case in a steady state situation 
even for very low Ca which should indicate that capillary 
forces dominate. 

If, on the other hand, Ca is large and viscous forces 
dominate, we will need a much larger difference in S to 
get a spanning cluster. This is because trapping of fluid is 
less likely to happen and therefore large clusters are more 
vulnerable to fragmentation through creation of bubbles. 
In the extreme case, where there are no capillary forces 
at all, a homogeneous mixing will occur. 



A. Pressure development 

For small Ca and subsequently large influence of capil- 
lary forces, the pressure build up is initially slow. A reg- 
ular structure of the fluid clusters will decrease the global 
pressure since the capillary part of the global pressure is 
dependent of the overall structure and the single phase 
pressure Ps Pc- 

For low capillary numbers the non-wetting fluid will be 
blocked in the area around the nodes due to the capillary 
barriers of the tubes. Wetting fluid will therefore in large 
parts of the system be immobilized. This trapping of 
fluid gives a lower fractal dimension D of the invading IP 
structure. 

For low Ca, the drainage of the IP part will eventually 
meet the compact imbibition front of the wetting fluid 
due to biperiodic boundary conditions. As long as the 
source of the drainage transient has an intact "belt" , as 
shown in Fig. 12 in the direction normal to the flow di- 
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Figure 3: Dynamic evolution of global pressure in a 128 x 128 
system at Ca = 3.2 x 10~^ and Af — 1. In the initial phase 
of the flow evolution the is a small build up in pressure, but 
when the the wetting front breaks through the initial belt of 
non wetting fluid the increase in pressure is violent. At the 
last, steady state, regime the fluctuations is approximately 
five times larger then in the first regime. 



rection, the invading front will move through bursts in 
the local pressure. These have a well defined distribution 
[l^. This regime is seen in Fig. pearly in the pressure 
evolution as a region with a slowly increasing mean pres- 
sure. 

When the first non-wetting fingers start to coalesce 
with the stable displacement of the wetting fluid after 
having completed a turn along the torus, a competition 
between two different processes arises. If the percolat- 
ing non-wetting cluster is not sustained and subsequently 
broken, Pc will increase. 

Eventually the initial IP cluster will break up and frag- 
ment into smaller clusters as shown in the third picture 
of Fig. 121 When there is no initial structure left in the 
system, the pressure saturates around a mean value with 
larger fluctuations than the previous case with the small 
bursts. 

The pressure development in Fig.|2|can be divided in 
three different segments. For this case where the injec- 
tion rate is slow, the flrst segment describes an IP regime 
where the Pc only describes local capillary fluctuations 
for the menisci movement along the front. The trapped 
clusters of wetting fluid will have little effect on the ef- 
fective viscosity of the non-wetting fluid and the front is 
not saturated. 

In the second segment we see a drastic increase in AP. 
Since the system is having the same injection rate and 
therefore is governed by capillary effects, this means that 
the effective dynamic viscosity finw has increased and the 
viscosity ratio M > 1. A qualitative explanation for 
this increase in viscosity can be seen in connection to 
the reduced permeability of the invaded, wetting region 
caused by trapped non-wetting clusters left behind by 
the compact wetting front. The non- wetting front has 
reached saturation width Ws- The result is a more stable 
front, and it was pointed out by Aker et al that Pc 
is proportional to the height difference in the front 

Pc oc Ah" , (8) 

where k — 1.0. 

In the third and last segment, the flow has reached 
steady state. The fluctuations in Pc are larger and the 
mean value higher due to the unstructured pattern of 
fluids. 



IV. STEADY STATE DISTRIBUTIONS 

In this section we look particularly at the distribution 
of non- wetting clusters. Because of the capillary barri- 
ers the non- wetting fluid has a much higher tendency to 
occupy pore-volume near the nodes rather than around 
the throats in the network. Due to volume conservation 
of wetting fluid, this will block the throats and create 
loop-less flngers of non-wetting fluid. Moving of trapped 
clusters must occur, and the steady-state fragmentation 
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and merging of non-wetting clusters N(s,t) is 



dtN{s, t) = [dtN] 

merging + [5t7V]ft.ag = . (9) 



A. Distribution of clusters 

The distribution of clusters depends on size. Fragmen- 
tation and merging of clusters in various dimensions are 
described i.e. in refs. [T ^ Iit I ITsI. 

At steady-state there will be a typical size s* of the 
clusters, and clusters of size s ^ s* are very rare. The 
distribution of clusters N{s, L) consists of both a regular 
and a singular part, N{s,L) = N,-eg(s,L) + 7Vsing(s,L), 
and for saturations Sn-w close to or equal to a critical 
saturation Sc, the singular part of iV(s, L) will dominate. 

Based on analogy to percolation and finite-size scaling 
|25l |. we assume the general expression for the singular 
distribution to be 



7V,ing(s,i)cxLXs-V(^) 



(10) 



where f{x) is a cutoff function decaying faster than any 
power law for a; ^ 1 and is constant for x <C 1. The 
exponent x which governs the system size dependence of 
iVsiiig, has the value x — d, in pure percolation. We have 
it here undetermined. 

When the clusters have a fractal dimension D, the typ- 
ical size s* oc where I = min(L, ^) where ^ is the cor- 
relation length. For large system sizes, the area A the 
clusters occupy is 



A = J sN{s,L)ds^ J s(iVrcg -f iVsi„g)ds 

— ^rog + ^sing ^ (H) 

and the singular part scales as 



/OO 
s^-^f{s/s*)ds 

L^{s*f-^ / x^-^f{x)Ax . (12) 



The lower limit of Eq. H12|) converges when we exclude 
arbitrary small clusters, and as a result 



when we use that s* oc near critical saturation. 

Long before steady-state, the non-wetting fluid forms 
IP patterns, and the total area A cx s*. The saturation 
of the invading, non-wetting fluid, then depends on the 
fractal dimension D and will scale locally as [l9| 



where L/oc is limited to the region where the defending 
fluid is purely wetting. Due to the incompressibility of 
the two phases, trapped fluids around the hull of the 
spanning, invading cluster causes the saturation of in- 
vading fluid to decrease even more as the system size L 
is increased. 

We have that A = Aging locally because this is in 
the pure IP regime, and there are no fragments of non- 
wetting fluid. Following Meakin the scaling of the 
singular cluster surface is 



A, 



siVsing(s,L)ds 



(15) 



which gives A^ing L''+(2-^)^ ^ and 



D + d 
D 



(16) 



In percolation, the fractal dimension of a spanning 
cluster is D{p = pc) = 91/48 giving r k, 2.05. In the 
regime where we have invasion percolation the area of 
wetting defending clusters behaves as 



A = j sN{s)A3 



(17) 



It was found that for IP structure, A grows as L'^, which 
imposed a r < 2 and A ~ ^ i^. This gives 

T = 1 + D/d and with a fractal dimension D = 1.82 
with trapped wetting clusters, this gives r = 1.91 for the 
distribution of wetting clusters. 

When steady-state is reached the mean saturation of 
the two phases stays constant over the entire system and 
every pattern created by the transients are wiped out. 
We write the non-wetting saturation as 



(18) 



where rig = L~'^N{s, L). According to Eq. (|ll|l we then 
get contributions to Snw both from a regular and singular 
part of Us- Since 5nw is constant and the spanning IP 
cluster is broken up, the probability that a given area- 
fraction of non-wetting fluid is attached to the spanning 
cluster changes from unity to 



Poo oc {Snw - ScY ^ L for Snw > Sc 



(19) 



D-d 
^loc ' 



(14) 



We take finite size scaling into account and the correla- 
tion length ^ cx \S — Sc\^'^ when Snw is chosen not to be 
far from Sc- 
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The singular contribution to Poo then gives the follow- 
ing relation between critical exponents 

!^ = {T~2)D + d~x>0 , (20) 
since /3 >0 and ly > 0, which gives 

r>2-^. (21) 

From this, we can have a value of r < 2, if the exponent 
X < d. In percolation where x = 2 we must have a r > 2 
in order to ensure a positive /3 > 0. 

When we look at the probability Us of having a cluster 
of size s in the limit of L ^ cx) it must approach an 
asymptotic, constant value. Since ris^sing < and 

rising - L^-'-s-^f (^) ^ Const. L^-^ , (22) 

as L — > oo, and hence s* — > oo, this means that x ^ 2 for 
finite s. It is clear the our situation differs from ordinary 
percolation. 

B. Evolution of wetting clusters 

It has been shown both experimentally and in numer- 
ical studies that the width of an invading non-wetting 
front depends both on Ca and the time t of the 
evolution. The front width w of the invading front is 
calculated as the standard deviation of number of points 
n{y) belonging to the front at a distance y from the av- 
erage position of the same front. It has been proposed 
that Wnw scales width time as 

w = t^''h{t, Ca) , (23) 

where h{t, Ca) is a crossover function. When the width 
of the front has reached saturation, i. e. at large time 
scales t ^ w^^^'', the saturation width Ws is no longer 
dependent of time and Ws is purely a function of Ca [2l| . 

At low Ca, the front is wide and trapped clusters of 
wetting fluid exists on many different length scales. This 
situation is best described as IP with trapping, with a 
fractal dimension D = 1.82 The distribution of 

these wetting clusters is shown in Fig. ^ We see a sig- 
nificant difference in the scaling at different times. The 
smallest exponent t according to Eq. H1U|) is r w 1.7 
which is consistent with the experimental findings by 
Frette et al At this point there is a cutoff for large 
clusters as these have not yet had the time to form. In 
the latter case r « 1.9 which is in more correspondence 
with the scaling proposed by Meakin [23| in Eq. (|17|l . 
When the system is in steady-state, the wetting clusters 
are no longer trapped and the largest wetting clusters 
will break apart. 
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Figure 4: Distribution of wetting clusters near the non- 
wetting front at different stages of evolution. The earliest 
stage (a) has a slope of 1.71. The cutoff here is due to lack 
of larger clusters that will evolve at later stages. As time in- 
creases (+), the slope is now 1.91, while at steady state (■) 
the power law distribution is completely overrun by the cutoff 
function. 

As Ca increases, the non-wetting front approaches 
a qualitatively more compact displacement |23l l24j . 
Looped fingers will also occur since the capillary barrier 
is weak compared to and there are no large, immobile 
wetting clusters as shown in Fig. |S| 

C. Effect of saturation 

For a certain saturation level of non-wetting fluid, 
at the equilibrium described in Eq. ^ the non-wetting 
fluid will percolate the system. We assume as in perco- 
lation theory the singular distribution of clusters to 
be 

iV,i„g(s)~s-^exp(-^) . (24) 

Since s* oc oc j^nw — Sc\^^^ , we can write the 
typical cluster size in means of saturation and get s* oc 
I'S'nw — 5c|~"^/'^. The exponent a = \/{vD), and for sim- 
plicity we define c cx 15* — Sc\^^'^ . 

We have performed numerical tests of the above pre- 
sumption with a matching viscosity ratio M = 1. The 
tests have been initiated with both a complete separation 
of wetting and non-wetting fiuids and a random mixing 
of the two phases. The question is whether or not the 
initial conditions will affect the steady-state outcome. 

If there exists a fluid configuration that stems from the 
combined process of drainage and imbibition, the out- 
come of the two initial configurations described above 
should have the same statistical properties. However, for 
low Ca we encounter difficulties connected to "locking" 
of immobile configurations due to large capillary barri- 
ers. The time-steps in our model become unphysically 





Figure 5: (Color online) Four different snapshots of configurations with different Ca taken at approximately the same time. 
In the upper left figure Ca = 3.2 x 10~^. For the next three figures Ca is much higher, Ca = 1, 10 and 100. The largest 
connected cluster is colored red in the online version of the article. The non-wetting front is now compact with high degree 
of fragmentation. The viscosity-ratio of the two fluids is M = 1, but the effective viscosity may be altered due to decreased 
permeability caused by fragmentation of the clusters. 



large. In order to compensate for these effects we apply 
certain shocks to the system where the input flow rate 
is increased drastically for a small amount of time steps 
compared to the total number of ran time steps. 

We now consider distributions of clusters obtained 
from different system sizes and initial configurations. 
The clusters are identified through a customized Hoshen- 
Kopelman algorithm [2^. The mass s of the clus- 
ters is a continuous variable, and in order to make a 



histogram we bin the identifed clusters with a binsize 
As — Sinax X 1/10000. The largest cluster Smax is five 
orders of magnitude larger than the unit size s = 1. 

As described in previous sections, the case where the 
two phases are completely separated gives intermediate 
configurations with different cluster distributions. The 
distributions depend on whether the fluid is invading or 
withdrawing from a region. In a study by Wagner et al 
[271] experiments and simulations of fragmentation of non- 



9 



10 - 

A lo' - 

z - 

10" r 
lo'r 



100 

s 



Figure 6: Cumulative distribution N{s' > s) of non-wetting 
clusters at steady state flow regime for 512 x 512 system at 
Ca = 3.2 X lO"'^. The results are from a simulation with 
Snw ~ 0.5 and the clear cutoff at high s indicates that the 
critical saturation Sc > 0.5. 



wetting fluid were carried out. The case here was to take 
a full IP pattern, and reverse the fluid current so that the 
flow regime was changed from drainage to imbibition. A 
result to note in Ref. is that when the non- wetting 
saturation of the IP pattern at breakthrough follows the 
relation in Eq. H14|l as expected, the flnal saturation of 
non-wetting fluid seemed to be independent of system 
size. 

In traditional simulation and experiments of imbibi- 
tion and drainage of different fluids, there is a inlet flow 
of one phase into an already fully saturated region of 
another phase. In our simulations, the effective satura- 
tion of the two phases are constant throughout the entire 
experiments. 

The percolation threshold for bond percolation in a 
regular square lattice is Pc = 1/2. The distribution of 
clusters in the system referred to in Fig.Elwith Snw = 0.5 
has a clear cutoff for large cluster sizes s. It is therefore 
clear that - in order to see clusters of non wetting fluid 
with size in order of the system size - we must have a 
larger saturation Snw and Sc > 0.5. 

With a non-wetting saturation Snw = 0.69 we see a 
tendency towards a power law distribution of cluster size 
s. A closer study at the results in Fig. [7| reveals a value 
of the exponent t — 1 = 0.91 ± 0.03 from the cumulative 
distribution N(s' > s) while a value of r = 1.93 ± 0.05 
slightly lower than 2 is drawn from the ordinary distribu- 
tion of N(s). This is not very far from the prediction of 
ordinary percolation and in consistence with predictions 
concerning similar dynamic models which indicate a 
T signiflcantly smaller than 2. 

It is hard to pinpoint the exact value of Sc in a system, 
but analysis of the cluster moments behaviour give us 
clear indications of the whereabouts of Sc- 
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Figure 7: Log-log plot of distribution of non-wetting clusters 
at Snw = 0.69. The upper figure shows the distribution N{s) 
for a system of size 1024 x 1024. The distribution N{s,L) 
exhibits a power-law with a slope of 1.93 ± 0.05. In the inset, 
a data-collapse of different system sizes is shown with a size- 
dependency of X = 1.8. The lower figure shows the cumulative 
distribution N{s' > s, L) for different system sizes Lx L going 
from L = 256 to L = 1024. The distributions have slopes 
around 0.91 ± 0.03. This combined with the data in the upper 
figure indicates a value of the critical exponent r = 1.92±0.04. 



D. Scaling analysis 

In order to study scaling behavior in the fragmentation 
and merging process of Eq. (Q, we examine different 
moments of cluster distribution. 

The fcth moment is defined as 



Mk = 



(25) 



and rig = ris^reg + '^s.sing, the moments will consist of 
both a r egul ar and a singular part. We use percolation 
methods |25l | to extract information of critical saturation 
and critical exponents. 
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Since the data we analyze is divided into discrete his- 
togram bins we substitute the integral of Eq. (|25|l with a 
sum and get for the singular part 

Mfc.sing = ^'^s (X L'^-^c^-i-fc , (26) 

s 

with the assumption that x — 2 < 0, where x was defined 
in Eq. ((TU|l. 

The total kth moment is ^3 

Mk^Ak + BkC+CW'-'', (27) 

where the last term stems from the singular part and 
dominates for fc > 2, since the finite-size contribution 
from Eq. 1)26(1 is small. 

In particular, the second moment M2 is dominated by 
M2^sing and 

= \S-S,\-^ , (28) 

when s* ^ near critical saturation. Then (x — 2)1/ 
gives a correction to the exponent 7 = (3 — t)/ct from 
ordinary percolation theory, and 7 = 7 + (x ~ 2);^. 

Results from different calculations of M2 for different 
values of Ca are shown in Fig. |H1 We notice that there 
are two regimes in the plot separated by a crossover at a 
certain value of jS'nw — ^'d- 

This can be understood finite-size scaling, where we 
near a critical density Pc apply the finite-size ansatz 




-0.2 -0.1 0.1 0.2 

S -5 „ 

nw f.Ca 

Figure 8: First moment of the largest cluster plotted as a 
function of saturations Snw for L = 512. The three curves are 
for hence Ca = 3.2 x 10"^, Ca = 3.2 x 10"* and ca = 3.2 x 
10~^. The percolation thresholds are shifted towards higher 
values of S„m as Ca increases. The transitions are located 
at approximately Sc = 0.675, Sc ~ 0.70 and Sc ~ 0.725 and 
hence collapsed. 



The critical saturations are moving towards larger val- 
ues as Ca increases. Since the fluctuations around the 
critical saturations are large, it is difficult to establish 
accurate value of Sc- However, Fig. El with the differ- 
ent values of S'cff obtained in Fig. |H1 gives a power law 
behaviour with 7 = 1.95 ± 0.06. This value is signif- 
icantly lower than the percolation value of 7 but this 
can be understood partly from the system size correction 
X-2^0. 



n = $[(p-p,)LV^] . (29) 
This implies that the effective, critical density PeS is 



Pcff = / P ( ^ ) dp 



(30) 



Hence, in our case we use this expression with our sat- 
uration for the variable p and get 



ScS — Sc 



A 



(31) 



Setting this new critical saturation in for Sc in Eq. 128|) 
we see the two regimes as 



My cx 



i\Sn 



Scl-" forl^n 



'I Const. for jS'nw 

when only Snw varies and L is constant 



Sc\ 3> 
Sc\ ■C pTu, 



(32) 




Figure 9: Log-log plot of the second moment for different Ca 
for L = 512. The data are collapsed using the values for Sc 
obtained previously. There is only a narrow range where the 
power-law is valid due to finite-size effects. For Snw < Sc far 
away from Sc we see a cut-off and for jS'nw — Sc\ -C l/L^^^ 
the function approaches a constant. However, the obtained 
value 7 = 1.95 ± 0.06 is reasonable. 
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Figure 10: Strength of largest non-wetting cluster for Ca ^ 
00. The data indicates a threshold at 5'nw = 0.91 which is 
lower than unity. We can explain this by a lower cutoff of 
small clusters which will be of significance when operating 
with no surface tension. The data are from systemsize L — 
512 and taken from the mean of ten different runs. 



Since the critical value of saturation of non-wetting 
fluid decreases with the capillary number, we expect more 
or less an approach towards percolation for extremely 
low Ca. On the other hand, when viscous forces are 
totally dominant and Ca — > oo, the expected threshold 
should be around unity. This expectation is based on 
the presumption that tiny bubbles of either phase will be 
overrepresented and effectively break up large clusters 
forming out of one fluid phase. Then no spanning non- 
wetting cluster can form even at very high Saw 

It is of course a matter of view how small we accept 
such bubbles to be, because at some point we must re- 
gard nodes to be connected even though microscopic frag- 
ments of wetting fluid separates them. These microscopic 
clusters, which we disregard in our calculations, do have 
an effect on Sc as can be seen in Fig ^| Here we see 
Sc = 0.91 which is significantly lower than unity. 

In order to find a value for cr, we use the ratio 



= i, for fc > 2 



(33) 



Then we can easily find an estimate for cr, and even 
though there is a finite-size correction to the scaling of 
Mfc, this will vanish when we measure the ratio in Eq. 

The measured value of a from the relation c oc |S'nw — 
Scl^^" is shown in Fig. [TTl We find a = 0.47 ± 0.08, but 
this value is dependent of the chosen value of the critical 
saturation Sc- This value is arguable. However, the ratio 
of 7 and a is assumed to be constant, and hence the 
value of T in our simulations. Since D — l/ai' we get the 
relation 



7(7 



D 



(34) 



Figure 11: Log-log plot of the ratio M3/M2 for L = 512. The 
solid line has a slope of —2.14 and gives 1/cr = 2.14 ± 0.08 
using Eq. This gives a = 0.47 ± 0.08. 



There is of course uncertanties connected to the size 
effects related to x~2, but with the data-collapse in Fig.[7| 
which suggested ax — 2^—0.2 and a fractal dimension 
for the non-wetting clusters in the area D = 1.82 — 1.89 
[2lll2^. the relation above should give a r « 1.97. This 
is consistent with our measured t = 1.92±0.04. 



V. CONCLUSION 

In this paper we have presented results of cluster dis- 
tribution in a two-phase flow network simulator in two 
dimensions. Different to most studies of porous flow sim- 
ulations in networks, we apply biperiodic boundary con- 
ditions. This allows us to study steady-state conditions 
which are more similar to those encountered deep inside 
natural reservoirs. 

The simulations show that there is a crossover from un- 
stable front propagation to a compact steady-state flow 
as the patterns formed by initial transients are broken 
up. From the pressure development in Fig. lathis change 
in flow regime will result in an increase in AP and a 
more violent evolution of burst avalanches. We believe 
that the fragments of non- wetting fluid sustain their frac- 
tal properties, but the Snw at which we get a spanning 
non-wetting clusters is dependent of Ca even though the 
fractal front is completely dissolved. 

The steady-state distribution of clusters is shown to 
share many characteristics with that of ordinary per- 
colation when it comes to critical properties and the 
extraction of critical exponents. We find a power-law 
distribution of non-wetting clusters near critical satura- 
tion. This is obtained regardless of the initial organi- 
zation of fluids outside the steady-state regime. Even 
though T = 1.92 ± 0.04 is slightly lower than the value 
for pure percolation it is consistent with results obtained 
from other dynamical models and fragmentation studies 
|30j |. From the order parameter Poo oc (5 — Sc)^, non- 
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negativity of the exponent /? > may be ensured even 
though T < 2. This is shown through a finite-size analy- 
sis. 
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